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ABSTRACT 


The  effective  v  iscosity  of  a  su^jcnsion  is  defined  to  be  the  four-tensor  which 


relates  the  average  deviatork  stress  to  the  average  rate  of  strain.  We  detennine  the 


effective  viscoaty  of  an  array  of  spheres  centered  on  the  points  of  a  periodic  lattice 

IS  tttte* . 


in  on  incompiessible  Newtonian  Buid^  The  formulation  involves  the  traction 
exerted  on  a  single  sphere  by  the  fluid,  and  an  integral  equation  for  this  traction  is 


derived.  For  lattices  with  cubic  symmetry  the  effective  viscosity  tensor  involves 
'^re- 

just  two  parameters.  -Tliey  are  computed  numerically  for  simple,  body-centered 


and  face-centered  cubic  lattices  of  ^hercs  with  solute  concentrations  up  to  9|l]|9^'or' 


I  •'  M 


iL'j.. u;- 

the  close-packing  concentration.  Asymptotic  results  for  high  concentrations  are 


. i 


obtained  for  arbitrary  lattice  geometries,  and  found  to  be  in  good  agreement  with 
the  numerical  results  for  cubic  lattices.  The  low  concentration  asymptdtk'etl 
sions  of  Zuzovsky  also  agree  well  with  the  numerical  results. 
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We  shall  calculate  the  effective  viscosity  of  a  suspension  of  rigid  spheres  in  an  incompressible 
viscous  fluid,  with  the  spheres  centered  at  the  points  of  a  periodic  lattice.  For  such  a  special 
configuration  we  can  obtain  results  for  all  values  of  thtiwiMMMaaKentration  c,  from  zero  up  to 
die  close-packing  value  c^. 

We  begin  by  proving  that  the  average  deviatoric  stress  in  the  suspenaon  is  linearly  related  to 
the  average  rate  of  strsdn.  The  effective  viscosiQ  fifju  is  defined  to  be  the  four-tensor  in  this  rela¬ 
tion.  It  is  determined  by  the  traction  of  the  fluid  on  a  single  sphere,  and  an  integral  equation  for 
this  traction  is  derived.  Its  kernel  is  a  periodic  Green’s  fvmction  which  can  be  represented  as  a  sum 
over  the  reciprocal  lattice.  The  integral  equation  also  involves  the  angular  velocity  of  the  qihere, 
for  which  a  second  equation  is  obtained. 

For  a  cubic  lattice  (simple,  body-centered  or  face-centered)  the  symmetries  of  the  kernel 
greatly  simplify  the  integral  equation.  It  can  be  shown  that  the  angular  velocity  of  each  sphere  is 
just  one  half  the  curl  of  r*--;  average  fluid  veloci^,  and  the  effective  viscosity  tensor  amplifies  to 

Myw  =  +  (1) 

Here  /i  is  the  viscosity  of  the  fluid,  Siju  is  one  if  all  the  subscripts  arc  equal  and  zero  otherwise, 
and  a  and  fi  ate  (Unctions  of  the  concentration  and  the  lattice  geometry  which  can  be  expressed  as 
integrals  of  the  traction  over  a  single  sphere.  We  use  a  Galerkin  method  to  solve  for  the  traction, 
following  the  procedure  used  by  2^ck  and  Honsy  (1982)  and  Zick  (1983).  The  results  for  a  and  fi 
arc  shown  in  Illustrations  1-6. 

For  low  concentrations  Zuzovsky  (1976)  and  Zuzovsky,  Adler  and  Brenner  (1983)  have 
obtained  asymptotic  formulas  for  a  and  ft  to  Our  numerical  results  corroborate  their  for 

mulas  and  determine  their  range  of  validity.  Frankel  and  Acrivos  (1%7)  obtained  the  leading  order 
term  in  the  asymptotic  expansion  for  a  at  high  concentrations,  but  only  for  a  simple  cubic  lattice. 
Zuzovsky  (1976)  corrected  their  result  and  also  gave  the  leading  order  term  for  ft.  We  have 
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obtained  additional  terms  in  these  expansions  for  a  simple  cubic  lattice  and  have  also  calculated  the 
corresponding  results  for  body-centered  and  face-centered  cubic  lattices.  Our  analysis  is  in  fact 
valid  for  arbitrary  lattice  geometries,  and  a  general  high  concentration  asymptotic  formula  for 
is  given.  The  asymptotic  formulas  for  the  three  cubic  lattices  are  in  good  agreement  with  our 
numerical  results  except  for  fi  in  the  face-centerd  case,  where  numerical  errors  at  high  concentra¬ 
tions  are  severe.  These  results  show  that  the  contention  by  Kiq)ral  and  Bedeaux  (1978)  of  a  singu¬ 
larity  at  less  than  close-packing  concentrations  is  incorrect. 

This  work  was  supported  in  part  by  the  Air  Force  Qff-trg  of  Sf.tentific  Research, 
the  Army  Research  Office,  the  nf  Waval  Research  and  the  National 

Science  Foundation.  Computing  time  was  provided  by  the  CLaSSlC  project 
of  the  Office  of  Naval  Research. 

Formulation 

We  consider  a  set  of  solid,  neutrally  buoyant  spheres  of  radius  b  centered  on  the  points  of  an 
infinite  three-dimensional  lattice  r®=aiai-l-oiai-f  ajaj,  aeZ’.  The  basis  vectors  at  deteimine  a 
unit  cell  having  volume  to  =  lai  faiXas)),  and  the  volume  concentration  of  ^heres  is 
c  =4wb^/3To.  We  assume  that  the  spheres  are  immersed  in  a  homogeneous,  isotropic  Newtonian 
fluid  undergoing  slow  flow.  Then  within  the  region  E  containing  the  fluid,  the  pressure  p,  viscosity 
p,  velocity  ui  and  stress  tensor  o/y  satisfy  the  Stokes  equations: 


=  0,  xi€E, 

dxi 


(2) 


^  duf  du, . 
c,j  =  -p8,j+p(-^+j^l  x,^E, 


da,j 


dx, 


=  0,  xi^E. 


We  introduce  an  overall  constant  shear  tensor  yy,  with  y«=0,  by  requiring  that  the  fluid 
velocity  i//  be  ihc  sum  of  a  linear  part  fijXj  and  a  part  which  is  lattice-periodic.  More  precisely, 
we  impose  the  condition 


>■  -.‘’3 

•  ;r 


•  .•I 


=  u,(x)  +  y(jrf.  x^E,  a€Zl 


(5) 


We  impose  no-sIip  boundary  conditions  on  the  surges  of  the  spheres.  The  spheres  may 
rotate  with  some  angular  velocity  Uf  which,  by  periodicity,  must  be  the  same  for  each  sphere.  Let 
5  denote  the  ^here  centered  at  the  origin  and  ni  the  outward  normal  to  its  surface.  We  can 
assume  that  5  is  not  transladiig,  hence  the  no-slip  boundary  condition  can  be  expressed  as 


=  ttjkUjXk,  jt*€9S’.  (6) 

Finally,  the  assumption  of  slow  flow  leads  to  the  vanishing  of  the  torque  on  each  sphere: 

I.„^tjkXjaynidA  =  0.  (7) 

oS 


Definition  and  Existence  of  the  Effective  Viscosity 


dU, 

Let  Stj  and  denote  the  volume  averages  of  the  stress  and  velocity  gradient,  respec¬ 
tively,  over  some  region  V.  As  shown  by  Batchelor  (1970),  they  are  related  by 


1  r  dU,  dUi  I  f 

■Stf  =  ^ J  [outXjm-n(u,nj+Ujn,)]dA.  (8) 

Here  K  is  the  volume  of  the  region  V,  the  sum  is  over  the  particles  in  V,  and  the  integral  in  die 
summand  is  over  foe  surfbce  of  the  summand  particle.  For  the  periodic  suspension  described  in 
the  last  section,  with  V  a  unit  cell  of  the  lattice  containing  S,  (8)  simplifies  to 


^iJ  =  +  To"‘/^0(*Xynt  dA,  (9) 

and  foe  bulk  deviatoric  stress  D/y  is  therefore  given  by 

^ij  =  +yji)  +  (ffa-vy  -  •»«)«*  dA.  (10) 
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Ftom  the  linearity  of  the  governing  equations  and  the  boundary  conditions,  the  integral  in 
(10)  nujst  be  linear  in  the  average  velocity  gradient  yij: 

J^autXjn^  dA  =  Mijuyu.  (11) 

Since  no  external  couple  is  acting  on  the  paiticles,  Mtju  must  be  symmetric  in  t  and  j,  and  for 
(tefiniteness  we  may  assume  that  =0.  rince  only  its  contraction  with  ym  is  signiScanL  By 
using  (11)  we  can  rewrite  (10)  in  the  form 


fl/y  =  2|i/y*/Ytf  02) 

where 

l^/JU  -  2t^(SitS^  +  SaSjk—jS(/Su)+JTfHAfljk/-J^tjJ^mmu)‘  (13) 

The  tensor  fi^  is  symmetric  and  traceless  in  its  first  pair  of  indices  and  traceless  in  its  last 
pair  of  indices.  In  Appendix  A  we  show  that 


I^IJU  =  FHtf*  ^^^) 

and  so  nyu  is  symmetric  in  k  and  /  as  well  Thus  the  average  deviatoric  stress  A;  is  linearly 
related  to  the  symmetric  part  of  yij,  which  is  die  average  rate  of  strain  tensor.  We  are  therefore 
justified  in  calling  nlju  effective  viscosity  of  the  suspension. 

The  formulation  becomes  somewhat  easier  to  use  if  we  multiply  both  sides  of  (13)  by  yu 
and  use  (11).  Then 


2p,}ii,Y/yT«  =  f*ir/y(Y/y+r/»)  +  To  *Y/y/  Xyai»n*d/i  (15) 

Now  the  effective  viscosity  niju  can  be  defined  as  the  unique  fourtensor  which  satisfies  (15)  for  all 
Y/y  with  Ya  =0,  and  which  also  satires  the  symmetry  conditions 
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a) 

b) M«u  =  0.  (16) 

c)  M/yu^Mw/- 

It  is  interestiiig  to  note  that  (IS)  is  precisely  the  foimula  one  would  get  using  the  energy  dissi¬ 
pation  method.  The  right  side  of  (IS)  is  the  actual  rate  of  cneigy  dissipation  within  a  lattice  cell, 
while  the  left  side  is  the  rate  of  energy  dissipation  that  would  occur  in  a  homogeneous  fluid  with 
viscosity  fbut^tensor  under  the  same  overall  shear  Y(/.  However,  the  symmetry  condition  (14) 
is  necessary  to  define  p^iju  uniquely  from  this  equation,  and  a  proof  of  this  identity  apparently 
requires  recourse  to  the  averaged  equation  formulation. 

Integral  Equation  Fomralation 

The  basic  unknown  quantity  in  equation  (IS)  ft>r  the  effective  viscosity  is  the  traction  // 
exerted  by  the  fluid  on  the  surface  of  the  sphere  5.  It  is  given  by 

//  =  X(^dS.  (17) 

In  this  section  we  use  Green’s  method  and  equations  (2)-(7)  to  derive  an  integral  equation  of  the 
first  kind  for  ft.  The  derivation  follows  that  of  Zick  and  Homsy  (1982). 

The  fiindamental  singular  solution  v^,  q^,  corresponding  to  unit  forces  in  the  ii:  direc¬ 

tion  ^plied  at  all  positions  y  +  r",  sqdsfies 


9  >»(»»?)  n 

dx, 

ruit(x.j)  =  -q*(xj)a(/  -I-  f*(■|“+■|^)’ 

OXj  ^ 


(18) 

(19) 

(20) 
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The  solution  for  is  (Hasimoto.  19S9) 


VitCxj)  = 


1 

4*^/iTo 


^Ik 

|S«|* 


g2iriS‘{x-j) 


Here 


S*  —  aibi  +  di2^  +  ttsbsi  <*  €Z^, 
is  a  vector  in  the  reciprocal  lattice,  which  has  basis 


-  hi  =  To'Hjhxas).  hi  =  To'HaaXai),  b,  =  xo'HaixajX 
The  prime  on  the  sum  in  (21)  means  that  the  index  a=0  is  to  be  omitted. 


(21) 


(22) 


(23) 


Let  r  be  a  unit  cell  containing  the  sphere  at  the  origin,  and  let  £o  =  be  the  fluid 
region  within  T.  For  y€£o,  we  multiply  (20)  by  U|(x)  and  integrate  over  £o  to  get 

/  M<(s)  dx  =  -«t(y).  y€£o.  (24) 

From  (2),  (3),  (18)  and  (19)  it  may  be  shown  diat 

dui  3vft 

Then,  by  using  (4),  (2S)  and  the  divergence  theorem,  we  can  write  (24)  in  the  form 
J^_^luirijt-v,i,iHj]njdA{x)  =  -i<»(y).  t 


(25) 


(26) 


We  now  proceed  to  simpUfy  this  equation.  If  /|  is  any  latrice-periodic  fonction  then  the  integral  of 

ftHf  over  the  surfoce  of  a  unit  cell  is  zero,  since  each  section  of  the  cell  boundary  has  a  congruent, 

periodkally  tran^ted  section  on  which  the  normal  points  in  the  opposite  direction.  The  tensors 

tj  dHMld  be  inieipreled  at  /  -/  ;  the  norawl  ib  poinu  out  oT  both  the  sphere  5  and  the  unit  cell  r. 

'*r-M  'ar  'm 


r^,  Vft  and  9^  are  all  periodic,  and  the  nonperiodic  part  of  ui  is  simply  yyXj:  Thus  the  integral 
over  dT  in  (26)  becomes 


=  /  7i,;t/Ty*nyd<4(x),  y€Eo.  (27) 

0T  97 

By  using  (19),  (20X  the  divergence  theorem,  the  periodicity  of  Vit,  and  the  fhct  that  yu  =  0,  we  can 
evaluate  this  latter  integral  The  result  is 

I^[uirijk-Vu,0(j\njdA(x)  = -yuyi,  y€Eo-  (28) 

The  integral  over  the  surface  of  the  sphere  in  (26)  can  also  be  simplified.  Using  (6).  (20)  and 
die  divergence  theorem,  we  can  write  die  first  part  of  it  as 

~  /  [«®«|Vjya~***itt/X*4(x-y)]dK,  y€E(^  (29) 

9s  S 

But  sO  since  ti^  is  symmeark  in  the  indices  i  and  ^  and  <(x-y)=0  because  ytS.  Thus 
this  integral  is  xero.  This  diet  and  (28)  enable  us  to  write  (26)  as 

+  = -M*(yX  y€£a  (30) 

We  now  let  y  tend  to  the  boundary  of  the  sphere  and  use  (6)  to  evaluate  the  fluid  velocity 
there.  Then  (30)  becomes 


X  »» wy fly  dA (x)  =  ytjyj  -  y€dS.  (31) 

This  is  an  integral  equation  for  the  traction  but  it  also  involves  the  angular  velocity  at. 

The  no*torquc  condition  (7)  provides  the  needed  supplementary  condition,  and  we  can  solve  (7) 
and  (31)  together  for  the  unknowns  /<  and  u,.  These  equations  are 


fvu(x,y)/t(x)dA(x)  =  y^jyj  -  eujUiyj.  y€3S, 

o5 

J^enkXjf(ix)dA(x)  =  a 


The  homogeneous  fonn  of  these  equations  has  the  solution  //(x)  =  kxi  for  any  constant  k. 
This  solution  is  everywhere  normal  to  the  surface  of  the  sphere  and  represents  a  uniform  pressure. 
This  presents  no  difficulties,  however,  for  since  yu=0  the  solvability  condition  for  (32)  is  satisfied, 
and  it  may  also  be  verified  that  any  uniform  pressure  in  the  traction  cancels  out  in  equation  (IS)  for 
the  effective  viscosi^. 

Symmetry  Relations  and  Cubic  Lattices 

For  simple,  body-centered  or  fece-centered  cubic  lattices  a  considerable  simplification  of  the 
problem  is  possible.  Their  sytrunetries  impart  a  structure  to  the  equations  which  allow  us  to  deter¬ 
mine  the  angular  velocity  u/  explicidy  and  to  reduce  to  two  the  number  of  independent  coefficients 
in  the  effective  viscosity  tensor.  To  show  this  we  define  //*  to  satisfy 

(34) 

as  ■’ 

and 

f  x,frdA{x)  =  0.  (35) 

95 

The  second  term  on  the  right  side  of  (34)  ensures  that  the  solvability  condition  for  the  integral 
equation  is  satisfied,  and  condition  (35)  pins  down  the  coefficient  of  the  null  solution  x/  in  //*. 
From  (32)  and  (34)  we  have,  up  to  a  solution  of  the  homogeneous  integral  equation,  that 

ft  =  (36) 

Consider  for  a  moment  an  orthogonal  linear  transformation  which  maps  the  reciprocal  lat¬ 
tice  onto  itself.  Then  fiom  (21)  we  can  show  that 

~  *’il»|('*pvyp)«  (37) 

and  by  changing  variables  from  x.  y  to  (x,  4y  in  (34)  and  (35)  we  can  subsequently  show  that 
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m 


mix)  = 


We  now  qsecialize  to  a  cubic  lattice.  With  no  loss  of  geneiali^  we  can  assume  that  the  com^ 
dinate  axes  and  the  principal  axes  of  the  lattice  coincide.  Then  tiw  reciprocal  lattice  is  also  cubic 
with  its  pimcipal  axes  coincident  with  the  coordinaie  axes.  Thus  it  is  invariant  under  the  reflection 
transfonnation  which  changes  die  sign  of  the  coordinate  of  its  argument,  and 

die  coortfinate  permutation  transformation  f  ^  =  iaO)j>  where  e  is  any  permutation  of  the  symbols 
{1,2,3}.  For  these  particular  otdK^nal  transformations  (38)  yields 


These  identities  can  be  used  to  simplify  (32)  and  (33).  We  first  solve  (33)  Ah'  W(.  Substitu¬ 


tion  of  (36)  into  (33)  yields 


q5 


However,  (39)  implies  that 


J^XykXjfTdA  =  i-lf  *****•'**“  J^ey^XjfrdA,  (42) 

and  since  dtis  must  hold  for  any  r  =  1,2,3,  the  integral  is  non-zero  only  if  the  indices  k,r,s  are  all 
diflierent  By  performing  a  coordinate  pomutation  «  taking  k-*!,  r-*2  and  s-*3,  and  using  (40), 
the  intcgrab  can  be  expressed  as 


S^XiJkXjfPdA  =  tknl^tijiXjfPdA. 


Upon  substitution  of  (43)  into  (41)  we  obtain 


iy n X rbUi)€ hrt  “0,  k  —  1,2,3. 
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The  soludon  of  this  equation  for  u/  is 

«/  =  -\iBjytj-  (45) 

This  is  just  one  half  the  curl  of  the  average  velocity,  which  is  the  same  angular  velocity  that 
would  be  exhibited  by  a  single  sphere  in  an  infinite  fluid  subjected  to  the  shear  yij.  Thus  in  the 
case  of  cubic  lattices,  each  sphere  spins  as  if  the  other  spheres  were  not  diere. 

The  structure  of  the  effective  viscosity  tensor  can  now  be  elucidated.  By  substituting  (45) 
into  (36)  we  obtain 

fi  =  il»»o(T/*+7jr)//*.  (46) 

which  upon  substinidon  into  (15)  yields 

^f^ljkiytjyu  =  PY<y(Yy +y^)  +  Tliyij(y„+y,r)LxjfrdA.  (47) 

qS 

This  can  be  simplified,  for  identities  (39)  and  (40)  readily  yield 

/  XjfrdA  =  Xj/TdA,  /  =  1,2,3  (48) 

qS  dS 

and 

I^^jfTdA  =  (49) 

Identity  (48)  implies  that  the  integral  is  zero  unless  the  indices  j,  i,  r,  5  are  equal  in  pairs,  and  (49) 
implies  that  many  of  the  remaining  non-zero  integrals  are  equal.  By  using  (48)  and  (49)  we  can 
write  relation  (47)  in  the  form 

2py«Y</Ytf  =  MY/y(Yy  +  Yy<)  +  YMYiy(Yu  +  +(6(y6*/  -  6/y«)/ H  (50) 

+(BikSji-8iju)lli  +i8iiSjk -8iju)[^], 

where 


The  solution  of  (SO)  for  niju  which  satisfies  the  symmetry  conditions  (16)  is 

where 


a  =  i(/tf-/il)  = 

(53) 

fi  =  jdu+ln)  =  jL  (xi/P  +X2/‘}^)dA. 

oS 

(54) 

Thus  the  effective  viscosity  tensor  involves  just  two  parameters  when  the  lattice  has  cubic 
symmetry.  They  can  be  found  from  the  solutions  to  two  vector  integral  equations  of  the  set  (34)- 
(35).  namely  those  for  //*  and  //“. 

Numerical  Method  for  Cubic  Lattices 

For  cubic  lattices  the  two  scalars  in  the  effective  viscosity  tensor  nJju  are  determined  by 
and  /“,  which  satisfy 

JiTo/  vy(x,y)//'(x)<//l(x)  =  yeOS.  r=l,2  (55) 

/  jt,/,i'(x)<//4(x)  =  a  (56) 

*dS 

These  equations  appear  to  require  a  numerical  solution.  The  major  difficulty  in  solving  (55)  is  the 
complexity  of  the  kernel  vy,  given  by  (21),  which  precludes  any  method  that  requires  numerous 
evaluations  of  v/y.  We  use  a  Galerkin  method,  with  the  basis  functions  chosen  to  be  analytically 
integrable  against  vjy.  This  method  yields  excellent  accuracy  with  relatively  few  lattice  sum  compu* 
tations  required. 
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Thc  Galerkin  method  may  be  outlined  as  follows.  We  write  in  the  form 

/i‘'(x)=  SarvH*).  (57) 

k,m 

where  is  a  set  of  basis  foncdons  on  dS  which  span  all  possible  solutions  for  The  basis 
functions  are  doubly  indexed  for  later  convenience.  Substitution  of  (S7)  into  (SS)  and  integration 
against  a  second  basis  function  yields  a  linear  system  of  infinite  rank  for  the  unknown  coefficients 

ar- 


i.'lAUgjar  =  bjr,  y  =  1.2.3;  /  =  1.2....;  /I  =  1.2,...  (58) 

/=!  kjH 

where 

=  f^f^li‘roVij(x,y)q>l"Hx)q>f''Hy)dA(x)dA(j)  (59) 

and 

V"  =  (60) 

The  expansion  (57)  is  truncated  and  the  resulting  finite  matrix  problem  is  solved,  yielding  an 
approximation  for  the  coefficients  in  (57). 

As  is  readily  verified  from  (21)  and  (55),  is  an  odd  function  of  x.  Thus  only  odd  func¬ 
tions  ^fl"^  need  be  included  in  the  set  of  ba^  functions.  Wc  choose  for  the  basis  functions  the 
polynomials 

ip^\x)  = 

9i^\x)  =  4*Jt2Xj*, 

S»i%)  =  x[*;t2*X3, 

where  {(pt,4t).  A:  =  1,2,...}  is  the  set  of  all  distinct  ordered  pairs  of  non-negative  even  integers. 
These  polynomials  can  be  derived  ffom  the  spherkal  harmonic  functions  and  form  a  complete  basis 
for  the  set  of  odd  functions  defined  on  the  surface  of  a  sphere. 


The  above  basis  functions  allow  the  integrals  in  (S9)  to  be  evaluated  analytically.  Since  the 


for  r+s  +  / =  iW+1,  0<r<s^  r  and  l^ic<  A/  (55^*'  is  also  needed  for  k=0.)  In  (65H67),  z 
denotes  lirbS".  For  A/ =  11,  which  corresponds  to  78  basis  functions,  there  are  only  684  such 
principal  lattice  sums. 

For  cubic  lattices  the  symmetries  of  v,j  and  the  basis  polynomials  produce  considerable  struc* 
hire  in  the  linear  system  (S8).  Specifically,  the  matrix  can  be  written  in  block  diagonal  form 
with  each  block  having  one  fourth  the  rank  of  the  full  matrix,  and  only  one  block  each  is  required 
to  c(»npute  the  coefficients  for  and  The  actual  systems  which  must  be  solved  to  compute 
and  //“  are 

^Uik  26^ 

aH  =  I  -bff  (68) 

A&k  A  fin 

and 

aSu  Aliu  aSu  oS  bff 

A}ht  Ajiti  A^  off  =  0  (69) 

Asiik  A^lt  Ayju  ® 

respectively.  Each  of  these  matrices  is  symmetric. 

The  matrix  in  (68)  is  singular,  with  null  vector  ai”=Si[i8i„  representing  a  uniform  pressure 
on  the  surface  of  the  sphere.  There  is  no  physical  basis  for  choosing  one  pressure  over  another 
(condition  (56)  is  merely  convenient),  but  as  previously  mentioned  a  uniform  pressure  does  not 
affect  the  effective  viscosity.  We  use  a  leas  squares  algorithm  to  obtain  a  solution  of  (68),  and  do 
not  implement  (56)  in  the  matrix  equations. 

The  solution  of  equations  (55H56)  is  not  our  final  goal.  It  is  rather  the  effective  viscosity, 
which  requires  the  integrals 


(xj/i^^dA  =  I  o^bB, 


f  =  1,2. 


(70) 


Thus  we  computed  the  effective  viscosity  by  solving  dte  linear  systems  (68)  and  (69),  conq>uting 


% 


i 


a  =  T2*S(a^i-aJD 


/J  =  T5:*S(aS+a8). 


and  substituting  these  into  the  formula  for  ityu>  repeated  here: 


y^'ju  =  /i(l+/S)*  j(5tt5y+5a5yt— -j5y6t/)  +  ^i(o-^X4y*/“yfiyfi*/X 


Asymptotic  Fonn  of  the  Effective  Viscosity  Tensor 


The  numerical  method  just  described  constituted  our  primary  method  of  solution  for  the 
effective  viscosi^  problem.  However,  as  will  be  illustrated  in  the  next  section,  the  numerical 
method  suffered  from  poor  convergence  properties  at  very  high  concentrations.  An  asymptotic 
analysis  for  the  high  concentration  situation  was  therefore  undertaken.  An  outline  of  the  analysis  is 
presented  in  Appendix  B;  the  results  are  reported  here. 

The  small  parameter  in  the  asymptotic  analysis  is 


e  =  l-(— )*^. 

Com 


The  final  result  for  the  equation  defining  the  effective  viscosity  tensor  is 

9  21 

Mtfwyyyw  ~  ItCmuYijyu  *"*^2  ninjnkm  (75) 

-  «/«/  “  +  <?(!)[ 

(I  h  h  ] 

Here  the  sums  are  over  all  unit  vectors  n  in  the  direction  of  a  nearest  neighbor  sphere,  with  respect 
to  a  single  reference  sphere,  and  the  matrix  Aij  is  defined  by 


-16- 


Aij  =  (76) 

» 

If  each  sphere  in  the  lattice  has  more  than  two  nearest  neighbor  H>heres  then  Afj  is  invertible.  If 
each  sphere  has  only  two  nearest  neighbor  qdieres,  in  directions  ±n,  then  A  is  singular  with  codi- 
mcnsion  1  and  null  vector  n.  For  such  two-nearest-neighbor  lattices  should  be  replaced  in 
(7S)  by  where  S  is  the  4x3  matrix  consisting  of  A  augmented  by  n^. 

When  (7S)  is  solved  for  n^u.  the  tenn  in  brackets  must  be  made  symmetric  and  traceless  in 
the  subscript  pairs  ij  and  k,l  (by  adding  terms  of  the  form  cijSu  or  BtjCu).  and  symmetric 
between  these  pain  (by  adding  a  term  of  the  form  where  Ciju  =  -cuij),  in  order  to  satisfy  the 
defining  conditions  (16).  The  precise  way  to  do  this  will  depend  upon  the  particular  lattice  being 
considered.  A  number  of  special  cases  have  been  examined  for  illustration,  and  the  resulting 
effective  viscosity  tenson  are  as  follows. 

L  Two-Nearest-Neighbor  Lattices,  with  incident  directions 

t^iju  -  0(1).  (77) 

II.  Cubic  Lattices: 

~  f»CmB[(^«"^-||-ln<)2(«<«y-7*tfXn*n/— jfiiw)  (78) 

~  "jlS  +  0(1) 

ni.  Tetragonal  Lattice  (orthogonal  basis  ,  with  |aii = |a2l  <  |a3|): 

Myw  -  fiCm«l(-^e-i-|^lne)2(«|fly- j«yX/i»«,- j««)  (79) 

~  ’^[BikBji  +  BiiBjii—j8ijB^-^SuBy+^8ij8uB^]lae]  +  0(1), 


where  Bij  =  '^ninj  ( n 


+ 


JJL  +JL 
N’“N 


). 


IV.  Hcxi^onal  Lattice:  Same  as  for  tetragonal  lattice,  but  with  the  fxtor  of  3/16  in  the  second 

line  of  (79)  replaced  by  1/8. 

Equation  (78)  for  cubic  lattices  can  be  written  in  the  form  (52),  from  which  asymptotic  for¬ 
mulas  for  a  and  fi  are  eaaly  obtained.  They  will  be  presented  in  the  next  section,  where  numerical 
estimates  of  the  leading  nonsingular  terms  in  the  asymptotic  expansions  will  also  be  given. 

Nnmeikal  Results  and  Discussion 

The  Galerkin  method  used  to  solve  (SS)-(S6)  involves  approximadng  by  a  finite  linear 
combination  of  basis  polynomials.  For  any  given  cakulation  all  basis  polynomials  of  degree  less 
than  or  equal  to  a  specified  M  were  used.  The  number  of  basis  functions  for  a  given  M  is 
(M  +  +2)/2,  but  because  not  all  basis  functions  contribute  to  each  the  xtual  number  of 

unknowns  is  2(M  +  IXA/  3)/8  for  /<*'  and  (A/  +  1X3A/ + 5)/8  for  /<“.  By  using  increasing  values 
of  M  while  keeping  other  parameters  of  the  problem  foced'^,  we  obtained  a  sequence  of  results  con¬ 
verging  to  the  exxt  values  of  a  and  fi.  Table  1  illustrates  the  convergence  of  the  solutions  as  a 
function  of  M  for  various  values  of  c.  These  results  are  typical. 

For  low  concentrations  convergence  is  very  rapid,  and  precise  results  can  be  obtained  with 
just  a  few  basis  functions.  Indeed,  for  concentrations  less  than  .06  only  3  basis  functions  (A/  :=  1) 
give  results  which  have  four  decimal  plxes  of  precision,  and  for  concenuations  less  than  .24  a  simi¬ 
tar  degree  of  precision  is  obtained  with  Just  21  basis  functions  (A/  =  5). 

For  higher  concentrations  the  convergence  rate  is  not  as  good.  As  c  increases,  the  traction  on 
tile  surface  of  a  sphere  becomes  more  and  more  peaked  near  points  which  arc  closest  to  other 
spheres,  and  more  basis  polynomials  are  needed  to  represent  it  accurately.  However,  the  matrices 
increase  in  size  with  the  number  of  basis  functions,  and  their  condition  numbers  increase  according¬ 
ly.  The  effect  of  the  inadequacy  of  the  basis  fuiKtions  in  representing  the  trxtion  is  magnified  due 

t  Since  the  matrix  for  ,1  given  M  rnnt.iins  all  the  clemcniN  nccctsarv  to  roim  the  matrix  for  a  smaller  M.  comput¬ 
ing  a  sequence  of  results  requires  only  slightly  more  ciTon  than  computing  one  result. 


to  these  high  condition  numbers.  This  is  particularly  noticeable  in  the  last  column  of  Table  1. 
where  the  computed  value  of  a  for  iV/  =  13  is  actually  less  than  that  for  M  =  11.  The  typical 
behavior  was  for  a  to  increase  with  M,  and  such  anomalies  were  never  seen  at  lower  concentra¬ 
tions. 

The  inadequacy  of  a  necessarily  finite  number  of  basis  functions  in  representing  the  traction 
at  high  concentrations  was  the  major  source  of  numerical  error.  There  was  strong  evidence  of  a 
downward  bias  in  the  high  concentration  computations,  for  all  three  cubic  lattices.  This  was  be¬ 
cause  the  basis  polynomials  could  represent  only  a  smoothed  approximation  to  the  traction  at  high 
concentrations,  and  so  the  computed  estimates  of  the  viscosity  coefficients  were  sometimes  too  low. 
The  face-centered  cubic  lattice  results  were  the  most  adversely  affected,  because  each  sphere  in  a 
face-centered  cubic  lattice  has  more  nearest  neighbors  than  spheres  in  cither  a  simple  or  body- 
centered  cubic  lattice.  Thus  the  traction  on  the  ^hcre  surface  is  more  complicated,  and  less  well 
representable  by  the  basis  polynomials,  for  face<entcred  cubic  lattices  than  in  the  other  two  cases. 

The  truncation  error  in  the  principal  lattice  sums  (6SH67),  os  measured  by  the  variation  in 
the  computed  coefficients  a  and  as  the  summation  range  was  increased,  was  insignificant  except 
for  high  concentrations  and  high  values  of  M.  Cancellation  and  roundoff  errors  in  the  pricipal  lat¬ 
tice  sums  were  deemed  insignificant  in  all  cases,  as  the  results  did  converge  with  increasing  summa¬ 
tion  range. 

The  speed  of  the  algorithm  depended  upon  the  type  of  lattice,  the  number  of  basis  polynomi¬ 
als.  and  the  summation  range.  The  computation  of  die  principal  lattice  sums  consumed  the 
preponderance  of  CPU  time,  and  the  complexity  of  this  computation  increased  as  the  cube  of  the 
summation  range  and  as  the  square  of  M  for  a  given  lattice.  The  sums  contained  twice  as  many 
points  for  a  body-centered  lattice  as  for  a  simple  cubic  lattice,  and  four  times  as  many  for  a  face- 
centered  lattice.  One  of  the  longer  runs  typically  required  several  minutes  of  CPU  time  on  a 
VAX-11/780,  while  a  low  concentration  run  required  only  a  few  CPU  seconds. 

As  previously  mentioned,  the  numerical  resuits  for  low  concentrations  arc  botlt  very  precise 
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and  quick  to  confute.  The  question  lemains  as  to  whether  the  results  arc  correct  Zuzovsky 
(1976)  obtained  low  concentration  asyiiqrtotic  results  for  what  we  have  identified  as  the  effective 
viscosity  tensor.  His  results  are 


o  -  ycll-(l-606)c-»-12ac*^-»-O(c’'^r‘. 


fi  ~  ■|c(l-(l-l-406)c-8ac*^-i-O(c’'^]-\ 


(80) 

(81) 


where  a  and  b  depend  on  the  lattice  geometry: 


SC 

BCC 

FCC 

a 

.2857 

-.0897 

-.0685 

b 

-.04655 

.01432 

.01271 

A  comparison  of  our  computed  results  with  (80)  and  (81)  is  illustrated  in  Table  2.  The  results 
agree  to  four  decimal  places  for  concentrations  of  .04  or  less,  with  tire  degree  of  agreement  slowly 
dropping  off  as  the  concentration  increases.  Zuzovsky’s  formulas  thus  confirm  that  our  program 
works  correctly  and  produces  accurate  results.  We  can  then  turn  the  tables  and  use  our  results  to 
check  his.  The  values  in  Table  2  and  similar  comparisons  for  the  other  lattices  indicate  that  (80) 
and  (81)  are  accurate  to  within  0.2%  for  concentrations  up  to  approximately  25%  that  of  close¬ 
packing,  and  to  within  S%  for  concentrations  up  to  approximately  50%  that  of  close-packing. 

For  high  concentrations,  we  ctunpare  the  numerical  results  with  our  asymptotic  results  for  a 
and  obtained  by  writing  (78)  in  the  form  (52).  These  formulas  contain  only  singular  terms  in  the 
small  parameter  e  =  l-(c/cmm)^^^  To  improve  them,  we  assume  that  the  leading  nonsingular 
terms  arc  1  and  elne.  Thus,  for  example,  we  write  the  difference  between  a  and  the  singular 
asymptotic  terms  for  a  simple  cubic  lattice  in  the  form 


a-(-^e-i-^lne)  =  C  +  D  elne  +  0{e), 
lo  oU 


(82) 


*.'.  V.  i:* A*.*^  «.■•«.".  c»  V’T*  ’rv  r-.  irj  T.' 
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where  C  and  2)  are  constants  to  be  determined.  By  plotting  the  computed  value  of  the  left  side  of 
(82)  against  cine,  we  determine  C  and  D  from  the  intercept  and  slope  of  a  line  drawn  through  the 
plotted  data  for  small  e.  We  do  the  same  for  all  other  a’s  and  fi’s. 

Unfortunately,  it  is  precisely  for  small  e  that  the  numerical  results  are  least  precise.  More¬ 
over,  in  computing  the  difference  between  the  numerical  and  asyniptotic  values,  we  may  lose 
significant  digits.  For  a  body-centered  or  face-centered  cubic  lattice,  the  graphs  for  a  exhibit 
definite  linear  behavior  as  e  becomes  small,  and  the  parameters  C  and  D  of  the  limiting  line  can 
be  determined  with  confidence.  But  in  the  other  cases  our  estimates  of  C  and  Z)  are  not  as  reli¬ 
able,  and  in  the  case  of  for  a  facc-centered  cubic  lattice  no  estimate  of  C  and  D  can  be  made. 

The  high  concentration  asymptotic  formulas  for  the  effective  viscosity  coefficients  a  and  /3, 
including  our  estimates  of  the  constant  terms  and  coefficients  of  e  In  e,  arc  as  follows: 

Simple  Cubic  Lattices: 

a  ~  -?7-e-*-^lne-l-3.1-025elne-l-O(e),  (83) 

iO  oU 

/J  ~  — J-  In  e  -I-  0.63  +  0.0e  Ine  -I-  0(e).  (84) 

Body-Centered  Cubic  Lattices: 

a  ~  In e  -  1.73  -  12.3 e  Ine  +  0(e).  (85) 

4 

/?-  -^e-‘-^^2.ine  +  i2.8  +  35elne-HO(e).  (86) 


Face-Centered  Cubic  Lattks: 


«  ~  •“«  +  5.7  -I- 15.5 e  In e  -»■  0(e).  (87) 

j2  ioU 

(88) 

Illustrations  1-6  show  our  computed  results  plotted  against  the  volume  concentration  of  the 
spheres,  and  also  show  graphs  of  Zuzovsky’s  low  concentration  asymptotic  formulas  and  our  high 
concentration  asymptotic  formulas. 

Appendix  A:  A  symmetry  relation  for 

In  this  appendix  we  prove  the  identity  where  nyu  is  the  tensor  appearing  in  (12). 

The  proof  is  a  modification  of  that  of  Hinch  (1972),  who  proved  a  related  result  for  a  single  particle 
in  a  fluid  undergoing  Stokes  flow. 

We  begin  with  the  reciprocal  theorem  of  Lorentz  (1896),  which  states  that  for  any  two  solu¬ 
tions  (u/.ojy)  and  (u),a)y)  of  the  Stokes  equations  (2)-(4)  in  a  volume  bounded  by  a  surface  F  with 
unit  normal  n, 

J^Uia’ij/tjdA  =  f^UfOtjnjdA.  (Al) 

Let  (ui,aij)  and  (u/.a'ij)  be  the  two  solutions  of  (2H7)  corresponding  to  y/j  and  y'lj,  respectively. 
Let  T  be  any  unit  lattice  cell  which  completely  encloses  the  sphere  at  the  origin.  Taking  F  to  be 
the  boundary  of  the  fluid  region  within  T,  the  first  integral  in  (Al)  becomes 

•Ct  dA  -  J^^eiuUkXia'ijnjdA.  (A2) 

Only  the  nonperiodic  part  of  uio'ij  contributes  to  the  integral  over  the  cell  boundary  in  (A2). 
Moreover,  the  no-torque  condition  (7)  implies  that  the  last  integral  in  (A2)  is  zero.  Thus  we  have 


•3' 

/I 

'  '^'  I 

m 

^.v*l 
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II 
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■>> 


:v>^ 
^;'j  ' 
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By  using  (3H6)  and  the  divergence  theorem  we  can  then  show  that 


J^_^u,a'tjnj  dA  =  iiToy,j(y,j+y'ji)  +  J^yfjXja'ucntdA. 
which,  using  (11),  may  be  written  as  . 


»to'ijHjdA  =  itToyij(y'ij+y’fi)  +  yijMijuyu. 

3T— 3» 


Similariy, 


/«  =  i^-roy’tjiyij+yjd  +  y'ijMijuyu- 

cT-^oS 


The  reciprocal  theorem  and  (13)  then  imply 


yjii^ijuyu  ~  yjti^ijuyu^ 


But  the  choices  of  yy  and  y'y  were  arbitrary,  up  to  tracelessness,  and  nyu  is  traceless  in  its  first  and 
last  piurs  of  indices.  Thus 


l^iju  =  PMty 


Appendix  B:  Asymptotic  Analysis 


We  consider  two  q)hercs,  each  of  radius  b  and  with  their  centers  a  distance  IL  apart  Let  n 
denote  a  unit  vector  along  the  line  connectng  their  centers.  We  define  a  set  of  cylindrical  coordi¬ 
nates  (r,B,z)  with  its  origin  at  the  midpoint  between  the  ^heres  and  its  ?-axis  in  tlte  n  direction. 
Let  ;fi  be  a  unit  vector  in  the  6=Q  direction  and  let  psfixm.  For  convenience  we  assume  that 

the  origin  is  translating  with  a  velocity  equal  to  the  average  velocity  of  the  spheres. 

* 

We  write  the  fluid  velocity  in  cylindrical  coordinates  as  u  =  i7f  4-  -i-  wz.  We  define  the 

parameter 


which  is  small  when  the  two  spheres  are  nearly  touching,  and  make  the  following  scalings: 


z  =  zLz, 
r  =  e'^Lr, 
u  =  Lu, 
y  =  Lv, 
w  =  t'^Lw, 

p  =  iu~^p. 

Equations  (2H4)  can  then  be  written  as 


z  =  t~^L~^z, 
r  = 

u  =  L~%  (B2) 

V  =  L"‘y. 
w  =  e~^L~^w, 

p  = 


Ur  +  r~^u+r~\+Wt  =  0,  (B3) 

-p,+u„  +  t[r~^UM-r-hr,-r-\-w„\  =  0,  (B4) 

-r“‘p,+ v„  +  t[r~^Ur0+r~^Ut+v„+r~\-r~^y-r~^W0,]  =  0,  (B5) 

-Pt  +  el  -  Mfi  -  r~\  -  +  r‘ V, + r"*»y„)  =  tt  (B6) 

The  boundary  condition  (S)  can  be  written  in  the  scaled  coordinates  as 


«±-je(/i^-l)«x  +  0(e^  =  ±(i>yA+piyijnj)an9  ±(-u'p+miyijHj)cos$  (B7) 

-  ezl-(«*/n)sintf  +  (u‘p)cos01, 

V ±-je(A^-  l)v, +(7(<^)  =  ±  (trp^mtfijnjysind  ±  (trih  +piyijnj}cos$  (B8) 

-  €^u-nr  +  ez[(wp)sin0  +  (u*>n)cos0], 

iy±-je(A^-l)w, +0(e*)  =  -(wrfi)rMntf  +  (u'^)rcos9  ±  e~'*ntyijnj  (B9) 

on  z  =  ±h,  where  A  =  l+r*/2. 


The  right  sides  of  the  boundary  conditions  can  be  divided  into  three  parts  corresponding  to 


the  relative  normal,  tangential  and  rotational  motion  of  the  spheres.  We  solve  (B3)-(B6)  for  each 
part  separately.  For  a  normal  motion  of  the  ^heres,  v  is  identically  zero  and  p,  u  and  iv  do  not 
depend  on  d.  We  express  these  latter  variables  as  3  regular  power  series  in  e,  and 

perform  a  standard  asymptotic  analysis.  The  resulting  differential  equations  can  be  integrated,  with 
the  boundary  conditions  determining  the  constants  of  integration,  until  a  last  second  order 
differential  equation  involving  the  pressure  is  obtained.  This  is  solved  with  the  following  regularity 
conditions:  (i)  p,  is  regular  at  r =0,  and  (ii)  p-*Q  as  r-*QO.  The  first  condition  is  necessary  so  that 
the  fluid  velocity  will  be  regular  at  /'=0,  and  the  second  condition  is  required  so  that  this  itmer 
solution  for  the  pressure  can  match  wifo  a  (regular  in  e)  outer  soludon.  A  similar  analysis  is  per^ 
formed  for  the  reladve  tangential  and  rotadonal  modon  problems  as  well.  The  combined  results 
are 


-  eUpiftjnj  +(m(ytjnj ~wp)cos0]h~^n  +  0{e^'^ 

u  =  + [(p/yy/iy+wm)sintf +(m;yy«y-wp)costf]A~^r 

+  e%iyijnj[{-^h~*r+3h~^rX2*-h*)  (Bll) 

+  (•^A~V--yA“V+-|-A~VXr"A^  -  ■l-rd-A*)]  +  0(e) 

V  =[i—miyijnj  +  u-p)sm6+(piyijnj+<a-m)cos$]h~^z-{u-n)rt'^  +  0{e)  (B12) 

w  =  e~\iyijnj[(h~^-^h~*)2^+jh~h]  (B13) 

+  [(piyijnj  +  u-m}^n6+{miyijnj~u-p)cos8}^[h~h^-l)r 
-  ((<o'm)sin^-(wp)cos^]r  +  0(e^ 


By  using  (3),  (17)  and  (B10XB13),  a  formula  for  the  traction  on  one  of  the  spheres  can  be 


fi  ~  (ib~^L  I  [e~^ntyijnji-^h~^r)  +  t~^l(ptyijnj  +  u)-ni)sind+(miyijnj-u-p)cos6]h~^ 

+  e-\y,jnji^h-^r+^h-'^r-^h-^r)  +  0(l)]r  (B14) 

+U~\yijnj^h-^  +  e~\yijnj(jh-^+^h-  +  jh-^ 

-t~'^(piyijnj+u-m)^nd+(miyijnj-u‘p)cos$]-2h~^r  +  0{l)]f 

The  dominant  contribution  to  the  integral  in  (IS)  at  high  concentrations  is  from  the  traction 
near  those  points  which  are  closest  to  neighboring  spheres.  By  multiplying  (B14)  by  Xj  and  in¬ 
tegrating  over  just  that  half  of  the  sphere  facing  its  neighbor,  and  summing  the  result  over  all  unit 
vectors  n  in  the  direction  of  a  nearest  neighbor  sphere,  we  obtain 

~  ^  •“«)  +  jStjntyun,  In e  (BIS) 

-2(5/*— /i/«*)y««/«y  Ine  -t-  leiuUkn/nj  Ine]  -f-  0(1). 

The  identity  ninj  +  mimj+PiPj  =  8ij  has  been  used  to  remove  any  explicit  dependence  of  (BIS)  on 
m  or  p. 

By  substituting  (BIS)  into  (7)  we  obtain  an  algebraic  equation  for  u*.  It  is 

=  ep;/Ytf  2  «;«/.  (B16) 

/} 

where  A  is  9ven  by  (76).  If  there  exist  only  two  nearest  neighbor  spheres  then  A  is  singular  and 
the  component  of  u*  along  the  line  through  their  centers  is  free.  However,  it  is  easily  verified  that 
such  a  free  component  docs  not  contribute  to  (BIS). 

The  final  equation  for  the  effective  viscosity  is  obtained  by  substituting  A^^tpjiya  2  tijiii  for 
u*  in  (BIS),  and  then  substituting  the  result  into  (IS).  After  reordering  the  subscripts,  and  using 
the  fact  that  To“*Z,’=3cniax/4ir,  we  obtain  (75). 
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Table  1:  Convergence  of  a  for  Simple  Cubic  Lattices 


.32 

.40 

.46 

.49 

2.6563 

4.3606 

5.8289 

6.5165 

2.8638 

5.3206 

8.4031 

10.479 

3.0151 

6.2231 

11.397 

15.353 

3.0239 

6.3882 

12.859 

19.388 

3.0256 

6.4440 

13.674 

22.260 

3.0257 

6.4541 

14.022 

24.499 

6.4541 

14.037 

24.256 

Table  2:  Numerical  and  Asymptotic  Results 
Low  Concentrations.  Simple  Cubic  Lattices 


c  1 

a  (comp.) 

a  (asymp.) 

B  (comp.) 

B  (asymp.) 

.012735 

.012735 

.012451 

.012450 

.01 

.025941 

.025943 

.024813 

.024812 

.02 

.053804 

.053810 

.049320 

.049316 

.04 

.11567 

.11570 

.97696 

.97677 

.08 

.26755 

26756 

.19337 

.19323 

.12 

.46580 

.46517 

.28995 

.28938 

.16 

.72502 

.72100 

.39009 

.38830 

.20 

1.0666 

1.0506 

.49665 

.49209 

.24 

1.5228 

1.4726 

.61306 

.60299 

.28 

2.1459 

2.0068 

.74379 

.72355 
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Figure  1 


Figure  2 
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Figure  6 


CAPTIONS  FOR  FIGURES 
Nunan  and  Keller 


a  for  a  simple  cubic  lattice  versus  concentration  c,  computed 
numerically  and  also  from  the  low  and  high  concentration  expansions. 

3  for  a  simple  cubic  lattice  versus  concentration  c,  computed 
numerically  and  also  from  the  low  and  high  concentration  expansions. 


a  for  a  body-centered  cubic  lattice  versus  concentration  c,  computed 
numerically  and  also  from  the  low  and  high  concentration  expansions. 


3  for  a  body-centered  cubic  lattice  verstis  concentration  c,  computed 
numerically  and  also  from  the  low  and  high  concentration  expansions. 


a  for  a  face-centered  cubic  lattice  versus  concentration  c,  computed 
numerically  and  also  from  the  low  and  high  concentration  expansions. 


3  for  a  face-centered  cubic  lattice  versus  concentration  c,  computed 
numerically  and  also  from  the  low  concentration  expansion. 
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